%     VIRTUAL FABRIC - программа для получения эффективных свойств текстильных композитов
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function [NODES_2D, ELEMENTS_2D, Min_Element_Edge, Max_Element_Edge, N_Layers_Suggestion] = Mesher_2D_Square(Points, Mesh_Seed_N, Element_Critical_Ratio)

global ErrCode;

global Merge_Tolerance_Multiplier;
global Auto_Thickness_Multiplier;
global Mesh_Seed_M_Multiplier;

NODES_2D = [];
ELEMENTS_2D = [];
Min_Element_Edge = [];
Max_Element_Edge = [];
N_Layers_Suggestion = [];

% Расчет числа элементов
Mesh_Seed_M = ceil(Mesh_Seed_N * Mesh_Seed_M_Multiplier);

Mesh_Seed_K = ceil(sqrt((norm(Points(1,:) - Points(5,:)) * norm(Points(3,:) - Points(6,:))) / ...
                        (norm(Points(10,:) - Points(6,:))*pi()/4 / Mesh_Seed_N * norm(Points(1,:) - Points(3,:)) / Mesh_Seed_N)));

% Проверка соотношения сторон элементов
Element_Ratio = norm(Points(1,:) - Points(5,:)) / Mesh_Seed_K / (norm(Points(10,:) - Points(6,:))*pi()/4 / Mesh_Seed_N);

if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка: недопустимое соотношение сторон элемента! [%f]\n', Element_Ratio);
    fprintf (1, 'Измените параметры построения сетки, либо увеличьте значение параметра Element_Critical_Ratio.\n');
    ErrCode = -1;
    return;
end

Element_Ratio = norm(Points(1,:) - Points(3,:)) / Mesh_Seed_N / (norm(Points(3,:) - Points(6,:)) / Mesh_Seed_K);

if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка: недопустимое соотношение сторон элемента! [%f]\n', Element_Ratio);
    fprintf (1, 'Измените параметры построения сетки, либо увеличьте значение параметра Element_Critical_Ratio.\n');
    ErrCode = -1;
    return;
end

% Начальные номера узлов и элементов
Next_Node_Number = 1;
Next_Element_Number = 1;

% Минимальная длина ребра элемента
Min_Element_Edge = min([norm(Points(9,:) - Points(8,:)) / Mesh_Seed_N, ...
                        norm(Points(9,:) - Points(10,:)) / Mesh_Seed_N, ...
                        norm(Points(6,:) - Points(9,:)) / Mesh_Seed_M, ...
                        norm(Points(3,:) - Points(6,:)) / Mesh_Seed_K]);
                    
% Максимальная длина ребра элемента
Max_Element_Edge = max([norm(Points(1,:) - Points(3,:)) / Mesh_Seed_N, ...
                        norm(Points(6,:) - Points(9,:)) / Mesh_Seed_M, ...
                        norm(Points(1,:) - Points(5,:)) / Mesh_Seed_K]);

% Построение сетки
[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Quadrilateral(Points(9,:), Points(8,:), Points(10,:), Points(7,:), Mesh_Seed_N, Mesh_Seed_N, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Quadrilateral(Points(15,:), Points(16,:), Points(14,:), Points(17,:), Mesh_Seed_N, Mesh_Seed_N, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Quadrilateral(Points(25,:), Points(21,:), Points(23,:), Points(20,:), Mesh_Seed_N, Mesh_Seed_N, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Quadrilateral(Points(29,:), Points(28,:), Points(33,:), Points(31,:), Mesh_Seed_N, Mesh_Seed_N, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

% -----

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(7,:), Points(8,:), Points(4,:), Points(5,:), Points(10,:), Mesh_Seed_N, Mesh_Seed_M, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(8,:), Points(9,:), Points(5,:), Points(6,:), Points(10,:), Mesh_Seed_N, Mesh_Seed_M, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(16,:), Points(15,:), Points(13,:), Points(12,:), Points(17,:), Mesh_Seed_N, Mesh_Seed_M, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(15,:), Points(14,:), Points(12,:), Points(11,:), Points(17,:), Mesh_Seed_N, Mesh_Seed_M, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(25,:), Points(23,:), Points(24,:), Points(22,:), Points(21,:), Mesh_Seed_N, Mesh_Seed_M, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(23,:), Points(20,:), Points(22,:), Points(19,:), Points(21,:), Mesh_Seed_N, Mesh_Seed_M, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(33,:), Points(31,:), Points(32,:), Points(30,:), Points(29,:), Mesh_Seed_N, Mesh_Seed_M, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(31,:), Points(28,:), Points(30,:), Points(27,:), Points(29,:), Mesh_Seed_N, Mesh_Seed_M, Next_Node_Number, Next_Element_Number, 1);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

% -----

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(3,:), Points(1,:), Points(6,:), Points(5,:), Points(10,:), Mesh_Seed_N, Mesh_Seed_K, Next_Node_Number, Next_Element_Number, 2);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(1,:), Points(2,:), Points(5,:), Points(4,:), Points(10,:), Mesh_Seed_N, Mesh_Seed_K, Next_Node_Number, Next_Element_Number, 2);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(2,:), Points(1,:), Points(11,:), Points(12,:), Points(17,:), Mesh_Seed_N, Mesh_Seed_K, Next_Node_Number, Next_Element_Number, 2);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(1,:), Points(18,:), Points(12,:), Points(13,:), Points(17,:), Mesh_Seed_N, Mesh_Seed_K, Next_Node_Number, Next_Element_Number, 2);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(18,:), Points(1,:), Points(19,:), Points(22,:), Points(21,:), Mesh_Seed_N, Mesh_Seed_K, Next_Node_Number, Next_Element_Number, 2);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(1,:), Points(26,:), Points(22,:), Points(24,:), Points(21,:), Mesh_Seed_N, Mesh_Seed_K, Next_Node_Number, Next_Element_Number, 2);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(26,:), Points(1,:), Points(27,:), Points(30,:), Points(29,:), Mesh_Seed_N, Mesh_Seed_K, Next_Node_Number, Next_Element_Number, 2);

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

[New_Nodes, New_Elements, Next_Node_Number, Next_Element_Number] = Mesh_2D_Segment(Points(1,:), Points(3,:), Points(30,:), Points(32,:), Points(29,:), Mesh_Seed_N, Mesh_Seed_K, Next_Node_Number, Next_Element_Number, 2); %#ok<ASGLU>

NODES_2D = [NODES_2D; New_Nodes];
ELEMENTS_2D = [ELEMENTS_2D; New_Elements];

% Допуск при слиянии узлов
Node_Merge_Tolerance = Min_Element_Edge * Merge_Tolerance_Multiplier;

% Слияние и перенумерация узлов
[NODES_2D, ELEMENTS_2D] = Merge_Nodes(NODES_2D, ELEMENTS_2D, Node_Merge_Tolerance);

% Рекомендация для числа элементов по толщине ячейки
N_Layers_Suggestion = ceil(norm(Points(10,:) - Points(29,:)) * Auto_Thickness_Multiplier / (norm(Points(10,:) - Points(29,:)) / (2*Mesh_Seed_N + 2*Mesh_Seed_M + 2*Mesh_Seed_K)));
